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We study the dynamics of particles in a multi-component 2d Lennard-Jones (LJ) flnid in the limiting 
case where all the particles are different (APD). The equilibrium properties of this APD system were 
studied in our earlier work [J. Chem. Phys. 142, 051104]. We use molecular dynamics simulations 
to investigate the statistical properties of particle trajectories in a temperatnre range covering both 
normal and supercooled fluid states. We calculate the mean-square displacement as well as angle, 
displacement and waiting time distributions, and compare the results with those for one-component 
LJ fluid. As temperature is lowered, the dynamics of the APD system becomes increasingly complex, 
as the intrinsic difference between the particles is amplified by neighborhood identity ordering and 
by the heterogeneous structure of the supercooled state. 


I. INTRODUCTION 

The tracking of individual particles and the study of 
their motions provides a simple means to probe the dy¬ 
namics of complex fluids, either experimentally or the¬ 
oretically (via computer simulations). The statistical 
analysis of particle trajectories can be used to uncover 
the physical origins of the broad distributions of relax¬ 
ation times and the heterogeneous dynamics (e.g., in par¬ 
ticle displacement distributions) observed in glass for¬ 
mers and colloidal suspensions close to glass and jam¬ 
ming transitions, respectively m, and in simple one- 
component fluids near the freezing transition in two di¬ 
mensions [S]. For example, computer simulation studies 
of glass-forming liquids suggest that (a) the slow dynam¬ 
ics and medium-range crystalline order are correlated 
(b) heterogeneous dynamics is a consequence of 
critical-like fluctuations of static structural order [5] , and 
(c) structural and dynamical heterogeneity are correlated 
m- Furthermore, particle tracking experiments on col¬ 
loids HIS] and theory/simulation studies of glass-formers 
and colloids (see e.g. ref. [5] and references therein) reveal 
that such complicated dynamical features are associated 
(at the single particle level) with intermittent particle 
motion- confinement of particle in cage formed by the 
neighboring particles and its subsequent release (trapped 
and jump-like motion); sting-like cooperative motion of 
particles is also observed mu]. 

Complex dynamics is also very common in crowded 
and heterogeneous bio-environments (e.g., living cells), 
as revealed by the single particle tracking experiments 
[TSHTT] . In living cells, many of the constituents (e.g., 
proteins) differ in shape, size and in their interactions 
with each other and with the surrounding medium; 
furthermore, they can be active (i.e., generate forces 
by converting chemical energy into mechanical work, 
upon ATP hydrolysis). There are several models of 
anomalous diffusion in biological systems, e.g., tran- 
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sient binding/unbinding, membrane compartmentaliza- 
tion, and membrane heterogeneity | 18 j (for a review on 
models of anomalous diffusion in crowded environments, 
see ref. [TO]). 

In the present study we consider a system of inter¬ 
mediate level of complexity between the former (glass- 
formers and colloids) and the latter systems (proteins in 
cells), and focus on the dynamics of particles with dif¬ 
ferent pair interaction strengths but with no size/shape 
polydispersity, that interact with each other and experi¬ 
ence random forces of thermal origin only. We will carry 
out a comprehensive statistical analysis of the dynamics 
of different types of particles in a multi-component 2d 
Lennard-Jones (LJ) fluid in the limiting case where all 
the particles are different (APD). Note that these parti¬ 
cles differ only in the strength of their attractive interac¬ 
tions with other particles and since these attractions be¬ 
come increasingly important as temperature is lowered, 
in this work we will focus on the liquid to solid transi¬ 
tion region, i.e. on the range of temperatures at which 
the differences between the particles are amplified and at 
the same time the system is not yet kinetically frozen. 
We will show that the dynamics of different particles in 
an APD system depends on their interaction parame¬ 
ters. In the fluid phase above the freezing transition tem¬ 
perature the motion of all the particles is diffusive and 
the only difference is in the magnitude of their diffusion 
coefficients; strongly interacting particles diffuse slower 
than weakly interacting ones. Qualitative differences be¬ 
tween the motions of different types of particles appear 
in the supercooled range below the transition temper¬ 
ature where the system becomes strongly heterogeneous 
(crystalline domains with mobile defects coexist with dis¬ 
ordered low-density fluid regions) and where the dynamic 
consequences of microphase separation between strongly 
and weakly interacting particles are revealed. 

Our paper is organized as follows: In section [llj we 
present the model and the simulation details and briefly 
describe the main results of our previous work on equilib¬ 
rium properties of APD systems. The analysis of mean 
square displacement, angle, displacement, and waiting 
time distributions is presented in sections |III| - |VI| In 
section IVTIl we summarize and discuss the results of this 






2 


work. 


II. MODEL AND SIMULATION DETAILS 


We perform MD simulations (in NVT-ensemble) in two 
dimensions with N = 2500 particles in a simulation box 
of size Lx = Ly = 60(t, where a is the particle diameter. 
At this number density p* = 0.6944 (= N/L^Ly) the 
system behaves as a liquid at sufficiently high temper¬ 
atures. Periodic boundary conditions are applied along 
both X and y directions. All the particles have same size 
a and mass m which are set to unity. Particles in the 
system interact via Lennard-Jones (LJ) potential 

Uij{r) = 4ey [(cr/r)i 2 _ (cr/r)^] , (2.1) 


where tij and r are the interaction strength and separa¬ 
tion between a pair of particles i and j respectively. The 
potential is truncated and shifted to zero at r = 2.5 (t. 
The equation of motion of particle i is given by the 
Langevin equation 


(Pvi dYi dU 

^ ~dYi 


fi , 


( 2 . 2 ) 


with Ti the position of particle i, C the friction coefficient 
(assumed to be the same for all particles), U is the sum 
of all the pair potentials Uij {i ^ j) for a given spatial 
configuration of the system, and fi a random external 
force with zero mean and second moment proportional 
to the poduct of temperature T and All the physi¬ 
cal quantities are expressed in LJ reduced units, with LJ 
time = 1 (the MD simulation time step used in the 
integration of the equations of motion is 0.005) [50] ■ The 
friction coefficient can written as C = with Td the 
characteristic viscous damping time which we fixed at 
50 and which determines the transition from inertial to 
overdamped motion (due to collisions with molecules of 
the implicit ’’solvent”) in the dilute system limit. At the 
rather high density in the present simulation, the damp¬ 
ing due to collisions between the particles dominates and 
the transition to the overdamped regime takes place on 
much faster time scales (of the order of 0.1). The simu¬ 
lations were carried out using LAMMPS code [21]. 

In order to characterize the interactions between the 
different particles, we have to introduce a distribution 
of pair interaction parameters that enter the LJ 
potential in Eq. |2.1[ In this work we follow the geo¬ 
metric mean (GM) prescription, introduced in ref. |22] 
according to which a particle i is assigned an interaction 
strength taken randomly from a uniform distribution 
in the range 1-4 (note that the identity of each particle 
is defined by its interaction strength and since the latter 
is a random number, all particles are different from 
each other). The pair interaction parameter between 
particles i and j is defined as the geometric mean of their 
interaction strengths , i.e., This leads to a 

peaked distribution of the pair interaction parameters 



FIG. 1: Distribution of pair interaction parameters P{eij) in 
the range 1 — 4 for GM system. 


with most probable value = 2.0 and mean 

(cij) = 2.42 (fig. [^. As a reference system, we simulate 
a one-component (1C) LJ system with the same density 
as the GM system and with pair interaction parameter 
tij = 2.5, the mid-point of the interval 1 — 4. 

The static properties of the GM system (as well as 
other probability distribution functions of pair interac¬ 
tion parameters) have been studied in detail in ref. [22] . 
The main finding of the above work was that as the 
system is cooled from the high temperature regime in 
which kinetic energy dominates over the attractive po¬ 
tential (and thus, the difference between the particles be¬ 
comes unimportant), the GM system undergoes a smooth 
transition into a self-organized fluid state characterized 
by short-range neighborhood identity ordering (NIO). In 
this state the identities of neighboring particles become 
strongly correlated in the sense that particles with high 
values of tend to cluster together (and similarly for par¬ 
ticles with low values of e^), resulting in microphase sep¬ 
aration between particles with high and low interaction 
strengths. NIO becomes increasingly more pronounced 
with decreasing temperature, a trend that continues un¬ 
til the system becomes kinetically frozen. 

It was found that the liquid-solid transition takes place 
within a narrow temperature interval in which the mean 
interaction energy changes rapidly upon cooling from the 
liquid to the “supercooled fluid” (using the terminology 
of ref. m) state, in which a hexagonally-ordered solid 
with mobile defects coexists with a rarehed fluid of par¬ 
ticles moving inside voids in the solid. The transition 
temperature can be defined as the temperature at which 
the radial distribution function decays over half the sys¬ 
tem size or by using other criteria. All methods yield 
T* « 1.0 for both the GM and the 1C systems (at the 
same density as in the present work). The precise nature 
of this transition is still controversial even for 1C systems 
and may be affected by finite size effects [531 [23] . Both 
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the static and the dynamic properties of the GM and 1C 
systems are expected to be quite similar for T ^ T* and 
differences between the two systems become significant 
as the temperature is decreased, approaching the upper 
limit (eij = 4) of pair interaction parameter values in 
the GM system, and neighborhood identity ordering be¬ 
comes increasingly important (there is no NIO in the 1C 
system). Since the differences between the two systems 
are amplified as the transition temperature (T* = 1) is 
approached, in the present study we focus on the temper¬ 
ature range around the transition temperature and com¬ 
pare the GM and 1C systems at the same value of the 
reduced temperature S = {T — T*)/T*. Note that the 
supercooled fluid regime extends from the equilibrium 
transition temperature {5 = 0) down to the temperature 
( (5 « — 0.5 for both the 1C and the GM systems) at 
which all the particles have condensed into the ordered 
solid phase and all voids and lattice defects became im¬ 
mobile on our simulation time scales. 

The longest time intervals in our simulations were 
2 X lO^r^j which corresponds to 40 million MD time 
steps. All the results reported in this work were ob¬ 
tained after the systems were equilibrated in the sense 
that ensemble-averaged properties (e.g., mean potential 
energy per particle) were time-independent and no fur¬ 
ther aging was observed. In order to get better statistics, 
all the distribution functions were calculated by collect¬ 
ing data along each trajectory and adding up data from 
the trajectories of all the particles in the particular en¬ 
semble studied. Note that averages obtained in this way 
involve both time and ensemble averaging and that the 
question of ergodicity (the difference between time and 
ensemble averages) did not arise. Nevertheless, we can 
not exclude the possibility of weak ergodicity breaking 
effects associated with probing individual particle trajec¬ 
tories over finite time intervals [5SJ [5^. The study of 
such effects will be the subject of future work. 

III. MEAN SQUARE DISPLACEMENT 

The simplest measure of the statistical properties of 
particle trajectories is the mean square displacement 
(MSD) during time t. We begin by calculating the MSD 
averaged over all the particles 

1 ^ 

(Ar2(t)) = — ^ (ri(t) - r,(0))^ , (3.3) 

i=l 

with i the particle index. 

In figure we compare the MSD for 1C and GM sys¬ 
tems (since all particles are different in the GM system, 
the average over all the particles means also that we aver¬ 
age over all particle types), at temperatures both above 
and below T*. Inspection of fig. [^a) shows that in agree¬ 
ment with expectation of simple diffusion in the normal 
fluid state, at temperatures above the transition d = 0.1, 
the MSD grows linearly with time for both systems. The 
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FIG. 2: MSD for 1C and GM (the latter is averaged over 
all particle types) systems shown at three different values of 
S indicated in the figure. Horizontal line shows (Ar^) = 1 
(particle size). 


temperature dependence of the diffusion coefficient in 
both systems for T > T* follows Arrhenius-like behavior 
and is shown in the SI. Note that the two MSD curves lie 
on top of each other above the transition temperature in¬ 
dicating that the corresponding diffusion coefficients are 
nearly identical. 

The difference between the MSDs in the 1C and the 
GM systems is resolved (the MSD for the 1C system is 
always lower than that for the GM system) at temper¬ 
atures below the transition, <5 = —0.1 and 6 = —0.2, 
where attractive interactions dominate over the kinetic 
energy and the effects of a broad distribution of interac¬ 
tion parameters survive the averaging over the particle 
types. Interestingly, as temperature is lowered through 
the supercooled range, the MSD curve in the 1C system 
develops a plateau at short times, followed by a transition 
to a subdiffusive regime {MSD oc at longer times. 
Upon some reflection we conclude that the plateau de¬ 
scribes the vibration of particles inside solid clusters with 
amplitude that corresponds to about 10% of the interpar¬ 
ticle separation (i.e., below the Lindemann criterion for 
melting), followed by escape from the “cage” at longer 
times, when the MSD reaches the particle diameter (the 
solid line in figs. (b) and (c)) and finally cross-over to 
subdiffusion, presumably due to the presence of mobile 
vacancies and defects. No short-time plateau is observed 
for the GM system, where the above-mentioned effects 
are smeared by the heterogeneous kinetics of different 
particle types. At very long times t > 10^, superdiffusive 
behavior is observed for the GM system at intermediate 
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temperature, S = —0.1 (see fig. ib)), possibly signify¬ 
ing a transition to yet another regime that lies outside 
the time scales reached in our simulation. The difference 
between the MSDs of the two systems disappears in the 
long time limit at yet lower temperature, S = —0.2 (see 

fig- He))- 



0 10 20 30 40 50 60 


FIG. 3: Typical spatial ordering of particles in GM system 
where the particles are colored according to their e; values, 
for (a) 5 = 0, (b) 5 = —0.1, and (c) 3 = —0.2. Color bar 
indicates the magnitude of Ci on a continuous scale from blue 
(1) to red (4). 


We now proceed to investigate the effects of particle 
identity (value of e^) on its dynamics in different temper¬ 
ature regimes. Typical configurations of the GM system 
at {6 = 0) and below (5 = —0.1 and —0.2) the liquid- 
solid transition are shown in figures (a) - [^c) where 
the color code indicates high (red) and low (blue) val¬ 
ues of Ci- Inspection of these figures shows that particles 
with larger values of Ci tend to lie inside solid clusters, 
while particles inside voids and at the boundaries and the 
defects of the solid are predominantly of smaller type. 

Neighborhood identity ordering is manifested by the 
magnitude of the effective interaction parameter of par¬ 
ticle i, defined as 

= («) 

with rib the number of nearest neighbor particles |22j. In 
fig. 0 we plot the ensemble average of this parameter, 
(ef^y^as a function of temperature, where the averaging 
is either over all the particles in the GM system (mid¬ 
dle panel), or over a given subset of all particles (upper 



FIG. 4: Variation of (e;®) with 3, shown for two different 
groups of particles with values of Si in the ranges 1 — 1.2 
(upper panel) and 3.8 — 4 (lower panel), and averaged over all 
particles in GM system (middle panel). 


and lower panels). The increasing splitting between the 
values of high and low particles shown in the up¬ 
per and the lower panels, respectively, is the signature of 
NIO: the tendency of particles with similar values of to 
cluster together leads to increase (decrease) of (e®®) for 
high (low) ei particles (at high temperatures NIO van¬ 
ishes and (e®®) —)■ y/2-5ei, in agreement with our simu¬ 
lation results). Surprisingly, fig. [^(middle panel) shows 
that NIO is a non-monotonic function of temperature 
below the liquid-solid transition, and decreases as tem¬ 
perature is lowered beyond 5 « —0.1. Inspection of figs, 
[^b) and (c) suggests this non-monotonic behavior is the 
consequence of condensation of (low e^) particles of the 
fluid on the solid clusters (composed of higher parti¬ 
cles) as <5 decreases from —0.1 to —0.2. Since the value of 
(e®®) depends on the local environment of the particle, 
this leads to the observed increase of (ef®) for the low 
Ci particles (lower panel in fig. and to the decrease 
of (e®®) for the high a particles (upper panel in fig. 
1^. The balance between the two effects is delicate but 
inspection of the middle panel in fig. shows that the 
latter effect dominates over the former one in the above 
temperature range. 

In order to distinguish between MSDs of particles with 
different values of e^, we consider four narrow (Ae^ = 0.2) 
intervals of Ci and plot the MSD as a function of time 
for these groups of particles (we average over trajecto¬ 
ries of particles in each group). As shown in figure [^a) , 
at the transition <5 = 0 (and at higher temperatures - 
not shown), the four groups of particles exhibit normal 
diffusion with diffusion coefficients that decrease with in¬ 
creasing interaction strength. If we define a relaxation 
time Tr as the time it takes to move one particle diame¬ 
ter on the average, we find that different types of parti- 
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time, t (in units of T^) 


FIG. 5: MSD for GM particles in four different ranges (as 
indicated in the upper panel), at temperatures (a) 5 = 0, 
and (b) S = —0.1. The intersection of the horizontal line at 
(Ar^) = 1 (particle size) with the MSD curves of particles 
with the largest and the smallest Si dehnes the range Att 
(see text). 


cles have relaxation times distributed over the range At^. 
with r,. ~ 1 for smallest particles and ~ 10 for the 
largest ones, as indicated in the figure. 

As temperature is decreased into the supercooled fluid 
range S « -0.1 (fig.[^b)), the dynamics of the GM sys¬ 
tem becomes heterogeneous in the sense that the simi¬ 
larity between the MSD curves of particles with different 
values of breaks down and At^ becomes much broader 
due to dramatic increase of the relaxation time for largest 
Ci particles (r^ ~ 100). In order to understand the origin 
of this dynamic heterogeneity recall that the transition 
temperature T* (« 1.0) reflects an average (over the dif¬ 
ferent particles in the system) property of the GM sys¬ 
tem; there are particles with > 2.5 (e^ < 2.5 ) for which 
the transition temperature of a system made of such par¬ 
ticles only, would be higher (lower) than that. For exam¬ 
ple, at r = T* a one-component system of particles with 
Ei —>■ 1 will not freeze because the transition temperature 
of such a system is about 0.4 (2.5 times lower that that of 
a IG system with = 2.5). Therefore, for T ^ T* one 
expects to observe crystalline clusters of larger parti¬ 
cles that coexist with a fluid of mobile smaller particles 
(see figs, [^a) and (b)). While this simple physical pic¬ 
ture is clearly oversimplified since there are small par¬ 
ticles that are trapped inside the solid clusters, it agrees 
with the observation of microphase separation between 
different types of particles in our system. The dynamical 
consequences of this segregation in the vicinity of T* are 
responsible for the fact that large particles diffuse rel¬ 
atively slowly compared to small ones, an effect that 


is clearly observed in fig. [^b) . 

We now focus on the dynamics of smallest particles 
at 5 « —0.1 (upper curve in fig.[^b)) which diffuse much 
faster and cover larger distances than other particles. 
Since the former particles are preferentially localized at 
the boundaries of the dense clusters, they can escape 
into the voids and move freely inside them. This free 
diffusion in a confined space may explain the observation 
of normal diffusion at short time scales followed by 
leveling out of the MSD due to the finite size of the 
voids at later times 10^ < t < 10'^. Beyond t « 10^ 
there is further increase of MSD that may be related 
to surface-tension driven coarsening of the voids on 
such long time scales. Glearly, the dynamics of the 
more mobile particles with the smallest values of is 
also responsible for the increase of particle-averaged 
MSD beyond t « 10"* observed in fig. [^b) (this mobile 
fraction gives the dominant contribution to large-scale 
displacements observed at very long times). If we now 
consider the dynamics of particles with the largest ei in 
fig. [^b) we notice that their MSD undergoes a transition 
from a plateau at short time scales (t < 10 — 100) to 
simple diffusion at longer times (t > 100), as they escape 
from a localized state characterized by lattice vibrations 
on length scales smaller than the particle diameter. For 
particles with intermediate values the behavior of 
MSD is in-between the two extremes discussed above. 

The above analysis gives us an crude picture of particle 
motions in our system. In order to obtain deeper insights 
about the dynamics of the APD system, we now pro¬ 
ceed to use other statistical methods of characterization 
of particle dynamics. One such method involves calculat¬ 
ing the distribution of relative angles between successive 
steps along the trajectory m- The sensitivity of this 
method has been demonstrated in the study of caging 
and escape in a quasi-2d colloidal suspension, near the 
jamming transition; in this regime the MSD curve is lin¬ 
ear, while the angular distribution shows a peak at 0 = tt 
( see fig. 5 in ref. m)- We proceed to calculate the angle 
distributions for 1C and GM systems at different tem¬ 
peratures near the solidification transition. 


IV. ANGLE DISTRIBUTION 


For a given particle trajectory, we define a vector 
V(<; A) = r{t -I- A) — r(t), where r(t) is the position 
of the particle at time t and A is the time lag. Then, 
the relative angle 9{t] A) between two consecutive vectors 
V(t; A) and V(t -|- A; A) along the trajectory is given by 
the dot product 


cos9{t; A) 


V(t;A).V(t + A;A) 
|V(t;A)||V(t + A;A)| ’ 


(4.5) 


where |V| is the magnitude of V. From eqn (4.5) we get 
9{t; A), and obtain the distribution of angles P{9; A) by 
computing all the relative angles along each trajectory 
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and averaging over the trajectories. A peak at 0 = 0 or 27r 
corresponds to correlated/directed steps, whereas a peak 
at 0 = TT represents anti-correlated steps, e.g. due to 
reflection from confining walls. By varying A we can ex¬ 
amine the dynamics at various levels of temporal coarse- 
graining. The upper bound for A is chosen such that 
A < t**, with t** the time to diffuse across the system 
size. In the following, we first examine the angle distri¬ 
butions obtained by considering the trajectories of all the 
particles, and then compare the angle distributions of dif¬ 
ferent types of particles by averaging over the trajectories 
of particles in a given interval Ae^ only. 



FIG. 6: Angle distributions P(0; A) of (a) 1C and (b) GM 
systems obtained by considering trajectories of all the parti¬ 
cles at (5 = —0.1. Different line types correspond to different 
values of temporal coarsening A as indicated in the inset. 


At high temperatures <5 > 0, the distributions are uni¬ 
form for both the 1C and the GM systems (not shown). 
As temperature is lowered below T*, one expects to ob¬ 
serve the effect of confinement by neighboring particles 
(caging) at short time scales, followed by escape from 
the cage and diffusion at a longer times. As shown in 
figure for 6 « —0.1 and small values of A, both 1C 
and GM systems exhibit a peak at 0 = tt in the angle 
distribution P(0; A) obtained by averaging over the tra¬ 
jectories of all the particles. With increasing A the height 
of the peak decreases and eventually the angle distribu¬ 
tions become isotropic (i.e., uniform in the range 0 — 27r), 
typical of normal diffusion where all the possible angles 
are equally probable. The difference in the dynamics of 
the two systems becomes evident when one examines the 
rate of change of P(0; A) upon increasing A; the decrease 
of peak height is much more gradual in the case of GM 
and even at A = 10^ the distribution is not completely 
flat. In the following, we will focus on the GM system 
and compare the behavior of P(0; A) for different particle 
types. 


< 

of 



e/2n Q/2n 


FIG. 7: Gomparison of angle distributions of the smallest 
and largest ti particles, shown for different A (top to bottom 
panel), at reduced temperatures S ~ —0.1 (left panels) and 
(5 « 0.1 (right panels). 


The angle distributions of two different groups of par¬ 
ticles with et in the ranges 1-1.2 and 3.8-4 respectively, 
are compared in figurefor three different values of A, at 
temperatures below and above the liquid-solid transition. 
At d « —0.1 and A = 10, a strong peak at 0 = tt is ob¬ 
served for large particles, as expected for particles con¬ 
fined within solid clusters. This peak decays slowly with 
increasing temporal coarsening and disappears (actually, 
it is replaced by a shallow dip) at A = 10"*^, indicative 
of slow cage relaxation. On the other hand, for small 
ei particles the confinement effect increases in strength 
at intermediate time scales and remains strong up to the 
largest time lag in our simulaton. Since such particles are 
expected to reside preferentially at the periphery of solid 
clusters and inside voids, this persistent behavior may 
be attributed to particles moving inside voids of different 
sizes which are forming and breaking as time goes on, 
with the largest holes being ~ 30cr in diameter (see fig. 
and supplementary figs. S2). Notice that from the MSD 
of the largest particles in fig. [^b), one might conclude 
that for t > 10 the particles are already well into the nor¬ 
mal diffusion regime, but this is not quite correct since 
the angle distribution shows that strong confinement ef¬ 
fects persist up to A ~ 500. This demonstrates that 
P(0; A) is a sensitive probe for dynamics which provides 
information that could not be obtained from the MSD 
alone. 


At higher temperature (S « 0.1) both particle types 
have uniform distributions in the range 0 — 27r; however, 
it is important to note that if we look very closely at 
the distributions for A = 10 and 10^, we realize that 
the curves do not exactly fall on top of each other. 
For instance, while at A = 10 the angle distribution 
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is completely flat for small particles, there is a dip 
at 0 = TT for large particles, indicating very weakly 
correlated motion. Such weak correlations can not be 
detected by analyzing the MSD curves, since all particle 
types have (Ar^(t)) ^ t when t ^ 10, even at 5 « 0 (see 
fig- da)). 

Although angle distributions provide information be¬ 
yond the MSD, they involve averaging over each trajec¬ 
tory (for fixed time intervals) and since the change of an¬ 
gle along the trajectory may be same for particles explor¬ 
ing large and small regions of space, some of the spatial 
information is lost. In order to provide detailed spatial 
and temporal information about the trajectories, we pro¬ 
ceed to study dynamic heterogeneities in the GM system 
using the displacement distribution. 


V. DISPLACEMENT DISTRIBUTION 


Distinct dynamical regimes are often observed in sys¬ 
tems close to criticality. The qualitative change in the 
dynamics across the transition is well captured by the 
displacement distribution (van Hove self-correlation func¬ 
tion) defined as 

Gsir, A) = {S{r - |r,(A) - ri(0)|)) , (5.6) 


where A is the time interva l, an d ri(A) is the position of 
particle at time A. Eqn (5.61 describes the probability 
of finding a particle at distance r from a given initial 
position after time interval A. 



FIG. 8: Comparison of the distribution of displacements at 
A = 500 for 1C and GM systems at four different values of 5 
indicated in the figure. 
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FIG. 9: Typical configurations of 1C and GM systems with 
particles colored according to magnitude of displacement r 
during time A = 500, shown for (a) 5 ~ —0.1 (upper panel), 
and (b) S « —0.2 (lower panel). Color code: gray open-circles 
for r < 0.75, green filled-circles for 0.75 < r < 1.5, and blue 
filled-circles for r > 1.5. 


In fig. we compare Gs(r, A) of 1C and GM systems 
at different reduced temperatures (5, for A = 500. The 
distributions for different values of A at fixed 5 (« —0.2) 
are shown in the SI fig. S3 At T > T* (S = 0.1) the 
displacement distribution is Gaussian (since the normal¬ 
ization is with respect to the measure dr, the Gaussian 
is multiplied by r), as expected for simple diffusion (see 
fig. I^a)). A qualitative difference between the the dis¬ 
placement distributions of two systems is observed at 
(5 « 0 (see fig. |^b)). Here, the displacement distribu¬ 
tion decays exponentially for 1C and is a complicated 
non-exponential function for GM, (see also supplemen¬ 
tary fig. [S4|(a)). For S <0, the distributions have several 
peaks: a primary peak in the range 0 < r < 0.75, and a 
secondary peak in the range 0.75 < r < 1.5, with the sec¬ 
ondary peak being more pronounced for 1C system (see 
figs.[8](c)-(d)). 

Insight into the origin of the multiple peaks of Gg 
for 6 < 0, can be obtained by coloring the particles 
according to the magnitude of displacement r. In fig. 
we use three colors that correspond to three different 
ranges of displacement: (1) r < 0.75 (primary peak), 
(2) 0.75 < r < 1.5 (secondary peak), and (3) r > 1.5, 
see figs. 1^ (c) and (d). Clearly, the first peak originates 
from vibrations of localized particles in solid clusters, 
the second peak arises from mobile particles associated 
with lattice defects and lastly, the higher order peaks 
that correspond to large displacements are due to mobile 
(fluid) particles located at the boundaries of the solid 
and inside voids. The relative change in the height of 
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FIG. 10: Displacement distributions of small and large ei par¬ 
ticles (left and right panels, respectively) are shown for three 
different values of A (see inset), at three different tempera¬ 
tures S (top to bottom panels). 


the peaks of Gg with changing d [compare figs. (c) and 
(d)] results from the corresponding change of the num¬ 
ber of particles in different r ranges in figs. |^(a) and (b). 


At longer times (larger A), the distributions broaden 
indicating that some particles are able to escape from 
their cages. The distribution decays much faster with r 
for large than for small ei particles. 


VI. WAITING-TIME DISTRIBUTIONS AND 
CAGE REMODELING TIME 


Inspection of particle trajectories at low temperatures 
shows that their motion is intermittent, i.e. consists 
of long periods of localized vibrations followed by sud¬ 
den jumps (see fig. II). Similar phenomena are ob- 


Tm 1281 [29] Inspection of the trajectories shown 
suggests that there is a broad distribution of 


served in a wide variety of systems that range from re¬ 
laxation in disordered materials to protein trajectories in 
live cells, 
in fig. [n 

trapping times (the amount of time particles spend in 
the cage, before hoping to another site). This leads to 
anomalous diffusion that can be described e.g., in the 
framework of the continuous time random walk (CTRW) 
model, in which a particle is trapped at a site for a ran¬ 
dom waiting time before jumping to the next site (for 
a discussion of this and other models of anomalous diffu¬ 
sion, see ref. [SO]). In this model, the probability density 
function (pdf) of follows the power-law form 


P{Tw) 


r-(l+a) 


(6.7) 


So far we discussed the displacement distributions 
obtained by considering all the particles in the sys¬ 
tem. In the following, we compare the distributions of 
different particle types in GM system. As before, we 
consider two different groups of particles with in the 
ranges 1 — 1.2 and 3.8 — 4 respectively, and calculate 
the corresponding displacement distributions (fig. 10). 
Above the transition {6 « 0.1) both particle types have 
qualitatively similar (Gaussian) distributions; however, 
the range is slightly larger for small particles. The 
distributions for small and large Ci particles are quali¬ 
tatively different at the transition (J « 0): Gs(r, A) is 
much narrower and decays more rapidly with r for large 
compared to small particles. In the latter case the 
form of the distribution is quite complex, especially at 
shorter times, suggesting that small particles sample a 
broader range of environments. It is interesting to note 
that the displacement distributions of small and large 
particles at <5 « 0, resemble the distributions obtained in 
a single-molecule tracking study of two different types 
of nuclear proteins Dendra2 and H2B, respectively (see 
figure 2 of figure supplement 1 in reference [HI). Indeed, 
the transcription factor Dendra2 behaves as a freely 
diffusing (and thus weakly interacting) particle, while 
H2B which is strongly bound to chromatin, displays a 
more restricted motion. At <5 « —O.I, the fraction of 
caged particles (with r < particle diameter) increases 
dramatically for all particle types and, as expected, the 
amplitude of the effect is largest for large particles. 


with 0 < a < I. Since the GTRW model has a distri¬ 
bution of waiting times and jump lengths, mapping an 
MD trajectory on a GTRW is non-trivial and requires 
filtering the continuous trajectories obtained from MD 
simulations |7]. However, in the present study we follow 
a simpler approach and calculate the time a particle re¬ 
mains confined within a certain threshold radius Rth [29j . 
Since depends on Rth one needs to check whether the 
exponent a is affected by the choice of Rth- In the fol¬ 
lowing we will take Rth = 1 which is about the mean 
inter-particle separation in the solid phase and thus cor¬ 
responds to the typical cage size. Therefore, P{tw) rep¬ 
resents the distribution of dwell times of particles in a 
cage (the distribut ion P{tw) for Rth = 5 is shown in the 
supplementary fig. |S5| ). 

For normal diffusion the waiting-time pdf corresponds 
to a Poisson process and is therefore exponential, 

PiTnj) = T~^ exp(-T„/r) , (6.8) 

with T the variance, see fig. [T^a). At <5 « 0, there is a 
qualitative difference between the waiting time distribu¬ 
tions of the two systems: while P{tw) of the IG system is 
still exponential, that of the GM system can be fitted by 
an exponential only at short times, < 10 (fig. [I2|[b)). 
Note that while in this range of waiting times, the pdf 
of the GM system decays faster with than that of the 
IG system, the situation is reversed at longer time scales. 
For (5 < 0, the waiting time distribution is well-fitted by 
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FIG. 11: Typical trajectories of particles in the bulk solid 
phase for 1C and GM systems at <5 « —0.2.The two GM 
trajectories correspond to particles with = 1.1 and 3.9, 
respectively. Particle positions are stored at intervals of 10 
and the trajectories are of length 2 x 10® (in units of rij). 



FIG. 12: Waiting time distributions P{tw) (obtained by con¬ 
sidering all particles) at three different values of 5, for i?th = 1- 
Exponential fits are shown by solid (r = 2.35) and dashed 
(t = 2.25) lines in (a) and by solid (r = 2.35) and dashed 
(r = 5) lines in (b), for GM and 1C systems, respectively. In 
fig. (c) the solid and dashed lines are power-law fits. 



FIG. 13: P(r„) compared for two different types of particles 
in GM system. Solid and dashed lines in (a) and (b) are 
exponential fits to the distributions for ti : 1-1.2 and ei : 3.8- 
4, respectively. The corresponding decay times are (a) 2.1 
(solid) and 2.8 (dashed), (b) 2.1 (solid) and 3.5 (dashed), (c) 
power-law fits. 



FIG. 14: Mean waiting time (r^,) as a function of 5 for 1C 
and GM systems. With increasing 5 (> 0) the value of (r^,) 
converges to the same value, irrespective of particle type. 


a power-law [eq. (6.7)] with a « 0.75 for 1C system and 
a « 1.3 for GM system (fig. [T^c)), indicating that in 
this temperature range escape from a cage is easier, on 
average, for particles in the GM system than in the IG 
system. 

An interesting question is whether different particle 
types in GM system have similar or distinct waiting time 
distributions. To this end we computed P{tw) for small¬ 
est and largest particles and compared them at differ¬ 


ent temperatures, see fig. 13 As one would expect, at 


higher temperatures (5 > 0, escape from the cage is al¬ 
ways a Poisson process and P(t^ ) is well-fitted by an ex¬ 
ponential distribution, eqn (6.8), though the relaxation 
time is shorter for low Ci particles. At the transition 
(5 « 0, the curves have a more complex form and a broad 
transition region from exponential to power law behavior 
(with a plateau in the range 10 < r^, < 30) is observed 
in fig. [I^b). At 5 < 0 both curves are well-fitted by a 
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13 


c). The fits yield 


power law P{tuj) ~ see fig. 

a « 0.75 for large and a « 1.8 for small Ci and the 
mean of these two values roughly corresponds to the ex¬ 
ponent a « 1.3 observed for GM system where we did not 
distinguish between particle types, see fig. [I^c). Com¬ 
parison of figs. fig. [^c) and[I^c) shows that in the low 
temperature range, the power-law exponent correspond¬ 
ing to high ei particles in the GM system is the same 
as that of the 1C system, suggesting that the dynamics 
of escape from a cage in a crystalline solid is dominated 
by (universal) excluded volume repulsions and does not 
depend on the strength of attractive interactions. 


A plot of the mean waiting time (tu,) as a function of 
the reduced temperature <5 is shown fig. Note that for 
(5 > 0 the value of (tu,) for the 1C system matches that 
of the particle-averaged GM system but for i5 < 0 the 
average waiting time for the 1C system coincides with 
that of the highest Ci GM particles. The above results 
concur with the previously made observation that the 
behavior of a particle in its effective cage is dominated 
by geometric (steric repulsion) effects that are similar 
for all the particles and depends only weakly on the at¬ 
tractive forces that distinguish between them. Thus, at 
high temperatures, the system is in a homogeneous liquid 
state and the effective cage is very similar for the 1C and 
the GM particles. Below the transition, most particles in 
the 1C system and the higher Ci particles in the GM sys¬ 
tem are in a crystalline (hexagonally ordered) state and 
again experience very similar “cages”. Below the transi¬ 
tion, GM particles with small Ci values find themselves 
inside and at the boundaries of voids and are therefore 
less constrained than in the liquid state, resulting in low, 
nearly temperature-independent values of (tw) . The ori¬ 
gin of the sharp peak of (tw) at 5 « 0 for GM system 
is unclear at present but it is clearly associated with the 
transition between the liquid and the solid cages that 
takes place at this temperature. 


Finally, we consider the distribution of cage remodel¬ 
ing times P{tr) where cage remodeling time G is defined 
as the time to completely change the identity of the cage 
formed by the nearest-neighbors of a reference particle 
(note that this time is much longer than the escape time 
from a cage). In figure [l^ a) we compare P(G) for 1C 
and GM systems, computed by considering all the par¬ 
ticles in the system. For d —>■ 0, the distribution P(G) 
decays faster for 1C system which thus has a smaller (tr), 
see fig.fl^a). There is abrupt increase in the value of (G) 
at the onset of freezing <5 —0+, see inset in figure [l^a) . 
Since particles in the GM system are different, the cage 
remodeling time depends on particle type: the pdf decays 
faster (slower) for particles with small (larger) e^, see fig¬ 
ure ll^b); consequently, as shown in the inset, the value 
of (G) is always higher for larger particles (at S « —0.1 
the difference between {tr) of smallest and largest par¬ 
ticles is more than an order of magnitude). Therefore, 
one can view the GM system below the transition as a 
system of mobile low particles diffusing inside the voids 
formed by the slowly rearranging solid clusters made of 



cage remodelling time, (in units of 


FIG. 15: Distribution of cage remodeling times P{tr) at two 
different values of <5 indicated in the figure: (a) for 1C and 
GM (all particles are included) systems, (b) P{tr) for highest 
and lowest ei particles. The corresponding mean remodeling 
times are shown in the insets (the mean remodeling time for 
1C is shown for comparison in the inset of (b)). 


larger a particles. 


VII. SUMMARY AND DISCUSSION 

We have studied the dynamics of particles in a multi- 
component fluid model where all the particles in the sys¬ 
tem are different in the sense that each particle has a dif¬ 
ferent interaction strength e^. The local self-organization 
of particles according to their identity becomes more 
pronounced as temperature is decreased throughout the 
fluid region, as shown in fig. |4 (see also reference [22)1. 
The only effect of this neighborhood identity ordering on 
the purely diffusive motion of individual particles above 
the liquid-solid transition temperature is that different 
particle types move with different diffusion coefficients. 
This dynamical “degeneracy” is lifted at and below the 
freezing transition as the coupling between long-range 
hexagonal order and local neighborhood identity order¬ 
ing (ordered solid domains are enriched by high parti¬ 
cles while low Ci particles proliferate in and around the 
low density regions (voids) and at mobile lattice defects) 
leads to heterogeneous dynamical behavior. 

Major differences between the dynamics of one- 
component systems and systems in which all particles are 
different are observed even when one considers ensemble- 
averaged quantities. At high temperatures, analysis of 
MSD shows that on average particle moves by simple dif- 
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fusion and that the diffusion coefficients are nearly iden¬ 
tical in one-component and APD systems (both exhibit 
Arrhenius-like behavior, with roughly the same activa¬ 
tion energy). Below the freezing temperature, the MSD 
curve of the GM system is consistently above that of the 
1C system, an effect that becomes more pronounced at 
late times due the presence of relatively large number of 
mobile particles inside voids in the GM system. Anal¬ 
ysis of angle distributions reveals the effects of particle 
caging in both systems, but the transition to uniform an¬ 
gle distribution (corresponding to normal diffusion) with 
increasing temporal coarse-graining is more gradual for 
GM systems. The displacement distributions of the two 
systems are qualitatively similar for (5 < 0 and <5 > 0, but 
are strikingly different at the transition (i5 « 0), with 
much slower decay of spatial correlations in GM system. 
The waiting time distributions are quite similar for (5 > 0 
and follow Poisson statistics, while for (5 < 0 both distri¬ 
butions have a power-law form but with different expo¬ 
nents: a « 0.75 (1C) and a ~ 1.3 (GM). 

The rich dynamics of the GM system is fully revealed 
when one groups together particles with similar interac¬ 
tion strengths (e^) and analyzes the trajectories of groups 
of particles with low and high separately. Indeed, as 
has been demonstrated in protein tracking experiments in 
living cells |16j , particles of different types moving in the 
same complex environment have very different dynami¬ 
cal histories, because of the different ways in which they 
interact with the surrounding medium. Thus, although 
particles of all types undergo normal diffusion in the liq¬ 
uid regime 5 > 0, their diffusion coefficients D depend on 
particle type: lower (higher) D for larger (smaller) par¬ 
ticles. This quantitative difference between the motions 
of different types of particles becomes qualitative at and 
below the transition, as revealed by the various meth¬ 
ods of analysis: MSD [fig. [^b)], angle [fig. [^, displace¬ 
ment [fig. and waiting-time distribution^fig. [Hb)- 
(c)]. Even though the way different particles move along 
their trajectories cannot be described in terms of sim¬ 
ple diffusion or sub-diffusion, using the above methods 
combined with snapshots of the particle configurations, 
we were able to construct a relatively simple physical 
picture of the dynamical phenomena that take place in 
supercooled APD fluids below the liquid-solid transition. 
Particles with higher interaction strength form large 
hexagonally-ordered clusters that contain mobile lattice 
defects and voids of different sizes; particles with lower 
interaction strength are preferentially localized in and 


at the boundaries of these voids and at lattice defects. 
While on short time scales high interaction strength par¬ 
ticles vibrate about their local equilibrium positions in 
the cage formed by their neighbors, at longer times they 
move around due to a combination of processes which 
include escape from the cage and rearrangement of the 
solid clusters due to the motion of defects and surface- 
tension driven coalescence of voids. Weakly interacting 
particles are enriched in the fluid phase at the boundaries 
of the solid and inside voids and therefore the “cages” in 
which they move reflect the distribution of void sizes. 
The separation of time scales allows us to propose a very 
crude but qualitatively correct picture of the heteroge¬ 
neous dynamics of the APD system in the supercooled 
regime: strongly interacting particles determine the ge¬ 
ometry of the space in which weakly interacting parti¬ 
cles move, with the result that the latter particles diffuse 
in a crowded and confined environment, reminiscent of 
geometry-controlled kinetics in ref. m- However, in our 
case this geometry changes slowly with time and it is 
this combination of different time scales that leads to the 
complex motions observed. 

At yet lower temperatures than those reported in 
the present study (around S « —0.5) both the one- 
component and the APD systems become kinetically 
frozen on simulation time scales, in the sense that all 
the particles have condensed into the crystalline phase 
and the defects are immobile. Note that since particles 
can no longer exchange their positions, NIO ordering is 
arrested as well and one obtains an “APD glass”. As we 
have shown lately in a study of the random bond model 
with particle exchange on a lattice |32j , this temperature 
is close to the onset of the annealed to quenched tran¬ 
sition below which different realizations of the quenched 
distribution of interaction parameters P(e-ij) yield differ¬ 
ent results. 
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SUPPLEMENTARY INEORMATION 

SI. TEMPERATURE DEPENDENCE OF DIFFUSION COEFFICIENT 



1/T 


FIG. SI: Diffusion coefficient D vs 1/T (Log-linear plot). Temperature dependence of the diffusion coefficient D{T). The inset 
shows the diffusion constant as a function of the rescaled temperature T\ = T/ < [T) >. 


The temperature dependence of the diffusion coefficient is studied for T > T*, see fig. |S1[ where the diffusion 
coefficient D is obtained from the ensemble average MSD. D(/T) shows Arrhenius-like behavior for large T: 

D{T) = , (SI) 

with Dq the diffusion coefficient at infinite temperature, Ea the activation energy, and ks the Boltzmann constant. 
For T ^ T*, D[T) deviates from the Arrhenius-like behavior indicating the onset of dynamics slowing down. From 
the fit in the Arrhenius regime we obtained Ea and find that both 1C and GM systems have roughly the same value 
of activation energy Ea « 1.95. Similar values of Ea for both the system is expected since the effective interaction 
parameter (e®®) which defines the local particle interaction is very close to that of the 1C system and weakly depend 
on r, see fig. (3) in reference 
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FIG. S2: Typical configurations of the 1C and GM systems compared at the same distance from the respective transition 
temperature. Shown in the figure for 5 = (T — — g.l, 0.0, and -0.1 (Top to bottom panel). 
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FIG. S3: Distribution of displacements Gs(r, A): The distributions for 1C and APD systems are compared at different 
values of lag-time A when S ~ —0.2. Inset is the zoom in showing secondary peaks. The presence of distinct dynamical regimes 
in the system is indicated by the secondary peaks. For small A the secondary peak is same for all the systems and it is enhanced 
upon increasing A. On average the particles in CM system cover relatively greater distance. 
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FIG. S4: Log-linear plot of the displacements distributions shown in fig. for (a) 5 « 0 and (b) S ~ —0.2. At 5 ~ 0, the decay 
of Gs is exponential for 1C, while it is slower for GM. For 5 « —0.2, both 1C and GM have the first peaks (corresponding to the 
vibration of particles around their mean positions) well described by the Gaussian, and the decay of the curves is exponential 
for 1C, while it is not quite exponential for GM. 
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6 = - 0.2 



FIG. S5: Waiting time distribution P{tw) at 5 = —0.2 (obtained by considering all the particles) is shown for i?th = 1 and 
5. Having different threshold radius do not affect the exponent a for 1C, while for GM the value of a is slightly lower for larger 
Rt\i- 










